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Thermal  Instability  and  Heat  Transfer  Characteristics 
in  Water/Ice  Systems 


YIN-CHAO  YEN 


INTRODUCTION 

The  phenomenon  of  convection  in  a  horizontal  water  layer  formed  by  melting  either  from 
below  or  from  above  is  of  special  interest  because  the  system  divides  itself  into  two  regions 
with  an  unstable  zone  lying  under  a  stable  layer.  This  is  caused  by  the  nonlinear  temperature- 
density  relations  of  water  possessing  a  maximum  density.  Figure  1  schematically  shows  the 
temperature  and  density  profiles  of  the  water  layer  continuously  formed  by  melting  ice  from 
below.  It  can  be  seen  from  curve  a  that  the  water  layer  is  stably  stratified  as  long  as  T ,  <  4°C 
and  it  is  unstable  when  T,  >  4°C  (the  upper  boundary  is  maintained  at  phase  transition  tem¬ 
perature,  i.e.  Ti  =  0°C).  It  is  distinctly  shown  that  the  water  layer  between  the  ice/water  in¬ 
terface  and  the  4°C  isotherm  is  always  stable.  However,  the  relative  height  (i.e.  the  ratio  of 
Z/H ,  where  z  is  the  distance  from  the  warm  plate  to  the  place  where  gmax  exists  and  //is  the 
total  layer  depth)  of  the  unstable  region  increases  as  T,  increases.  Figure  2  is  a  simple  inver¬ 
sion  of  Figure  1  and  shows  the  case  of  ice  melting  from  above.  In  this  case  T,  is  at  the  melting 
temperature  (i.e.  T,  =  0°C)  and  Tt  >  0°C.  In  contrast  to  Figure  1,  the  relative  height  of  the 
unstable  layer  decreases  with  increasing  Tt.  Therefore,  for  water  with  its  maximum  density 
near  4°C,  even  for  a  constant  temperature  gradient,  the  density  profile  is  a  nonlinear  func¬ 
tion  of  the  vertical  ordinate  z. 

These  situations  are  quite  different  from  the  classical  Rayleigh  problem  in  which  a  con¬ 
stant  temperature  gradient  and  a  linear  temperature-density  relation  of  the  fluid  are  assumed. 
In  the  classical  problem,  a  fluid  layer  of  constant  depth  confined  rigidly  from  top  and  bot¬ 
tom  but  of  infinite  horizontal  dimensions  is  either  cooled  from  above  or  heated  from  below, 
and  convection  occurs  if  the  temperature  gradient  exceeds  a  certain  critical  value.  The  criter¬ 
ion  is  usually  expressed  in  terms  of  a  Rayleigh  number  with  a  critical  value  of  =  1708,  which 
has  been  verified  by  Blnard  (1900)  who  observed  the  regular  cell  structure  for  Ra  >  1708. 
This  classical  problem  has  been  studied  intensively  since  the  early  works  of  Rayleigh  (1916) 
and  Jeffreys  (1926).  Chandrasekhar  (1961)  provides  a  comprehensive  analysis  and  summary 
of  these  studies.  The  aim  of  this  study  is  to  review  and  analyze  the  stability  problem  due  to 
the  density  anomaly  of  water  and  its  effect  on  the  heat  transfer  characteristics  in  ice/ 
water  systems  due  to  the  evolution  of  the  interfacial  velocity  that  results  from  the  density  dif¬ 
ference  between  ice  and  water. 


ANALYTICAL  STUDIES  ON  THE  ONSET  OF  CONVECTION 
IN  A  HORIZONTAL  WATER  LAYER 

The  influence  of  the  density  anomaly  of  water  on  the  onset  of  convection  has  been  report¬ 
ed  by  quite  a  few  investigators.  Veronis  (1963)  was  first  to  investigate  the  phenomenon  of 
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Figure  1.  Schematic  representation  of  the  density  profile  in 
the  conduction  state  for  upper  boundary  T2  =  0°C  (equiva¬ 
lent  to  melting  from  below). 


Figure  2.  Schematic  representation  of  the  density  profile  in 
the  conduction  state  for  lower  boundary  T,  =  0°C  (equiva¬ 
lent  to  melting  from  above). 


penetrative  convection.  By  representing  the  water  density  near  4°C  with  the  quadratic  ex¬ 
pression 

Q  =  ^mau)2!  (I) 

where  emax  is  water  density  corresponding  to  Tmax(-  3.97°C)and7  -  7.68  x  10'*  and 

using  linear  stability  theory,  Veronis  summarized  the  computed  critical  Rayleigh  numbers 
Rac  for  both  free  (Table  1)  and  rigid  (Table  2)  boundaries.  In  both  tables,  X  is  the  ratio  of  the 
total  layer  depth  H  to  the  4°C  layer  depth  d  in  the  conductive  state  (i.e.  X  =  H/d),  a  is  the 
wave  number,  and  Rac  is  the  critical  Rayleigh  number.  In  Table  2,  n  is  defined  as  the  ratio  of 
the  rotation  rates  of  the  outer  to  inner  cylinder  [(1  ~n)  =  X],  Tac  is  the  critical  Taylor  num¬ 
ber.  The  values  of  fl!/X'  can  be  compared  with  values  of  c!/(l-^)!7r:  and  Rac/2X47r4  with  Tac/ 
2(l-u)4ir4.  There  is  of  course  a  quantitative  difference  between  the  results  of  the  two  cases  be- 


Table  1.  Critical  Rayleigh  numbers  and 
corresponding  wave  numbers  for  vari¬ 
ous  values  of  X  in  rigid-free  boundary 
problem  (after  Veronis  1963). 


Table  2.  Critical  Rayleigh  numbers  and  corresponding 
wave  numbers  for  various  values  of  X  in  rigid-rigid 
boundaries  (after  Veronis  1963). 


cause  of  the  different  boundary  conditions,  but  it  can  be  observed  that  the  qualitative  behav¬ 
ior  is  the  same  in  both  cases.  In  both  tables,  the  unstable  layer  height  d  is  used  in  defining  the 
Rayleigh  number.  For  the  specific  cases  of  T2  =  4°C  and  8  °C,  i.e.  X  =  1  and  2,  Veronis  re¬ 
ported  Rac  values  of  13.437T*  and  87.187T*  respectively.  These  tabulated  values  of  Rac  clearly 
demonstrate  that  the  onset  of  convection  is  dependent  on  the  thermal  boundary  conditions 
imposed. 

Sun  et  al.  (1969)  extended  this  work  by  adding  a  cubic  term  to  the  temperature-density  re¬ 
lation 

Q  =  emax[1-7i(7'-7max)2-7!(7'-  ^max)’]  (2) 

that  is  valid  up  to  30°C.  Following  a  procedure  commonly  employed  in  linear  stability  analy¬ 
sis,  a  modified  Rayleigh  number, 


2-ylAg(A7)’H'(l+  y^  A  AT) 


is  derived  in  which  A  =  (T>-Tm3X)/(T-T2),  H  is  the  total  melt  layer  depth,  and  AT  is  the 
overall  temperature  difference  (i.e.  AT  =  T,-T2).  They  found  the  Ra  values  were  functions 
of  two  parameters  defined  as  follows: 
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Figure  3.  (Rac)Xt=0  as  a  function  of  X, 
(after  Sun  et  al.  1969). 
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Figure  4.  -  1  vs  with  X,  as  a 

'Kac)\,=0 

parameter  for  the  rigid-free  case  (after 
Sun  et  al.  1969). 


Rac  values  are  computed  for  ranges  -4.25  < 
X,  <  -0.5  and  -1.4  <  X2  <  1.6.  Numerical 
values  of  (RaJX;  =  0  (i.e.  for  the  case  of 
parabolic  temperature-density  relation,  72 
=  0,  and  X.  =  -l /A  as  a  function  of  X,  were 
reported  and  shown  in  Figure  3  for  both 
rigid-free  and  rigid-rigid  boundaries.  The 
Rac  values  for  cubic  and  parabolic  tempera¬ 
ture-density  relations  were  expressed  in  a 
ratio  as  (RaJ^/fRaJ^  =  0  and  was  shown 
vs  X2  using  X,  as  the  parameter  for  rigid-free 
(Fig.  4)  and  rigid-rigid  boundary  conditions 
(Fig.  5).  These  results  are  valid  as  along  as 
the  melt  layer  contains  a  density  extremum 
within  the  boundaries. 

Merker  et  al.  (1979)  studied  this  same 
problem  with  the  representation  of  the  tem¬ 
perature-density  of  water  approximated  by 
three  different  polynomials  having  2,  3,  and 
5  terms.  Linear  stability  analysis  was  used 
and  the  resulting  perturbation  equations 
were  solved  with  the  aid  of  Galerkin’s  meth¬ 
od.  Figure  6  shows  the  general  diagram  of 
stability  for  fluids  having  a  density  extrem¬ 
um  in  which  the  values  of  d,  and  are  the 
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ter  Sun  et  al.  1969). 
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Figure  7.  Critical  Rayleigh  number  Rac 
as  a  function  of  T,  with  T  2  as  a  param¬ 
eter  for  Tw  (wall  temperature)  -  con¬ 
stant  (after  Merker  et  al.  1979). 


|°5  r  i  m  i  |  \  |  1  I  i  |  i  3  coefficients  of  thermal  expansion  evaluated 

r'y  \  q.=const  -  at  y  an£j  respectively,  and  the  values  of 

~  /  \\  „,5  -  Ra,  and  Ra2  are  the  corresponding  Rayleigh 

/  8/  \°  \t2  =  4°c  ~  numbers,  defined  as  Ra  =  gH'!3AT/av  in 

io“  —  \  \  — ^  which  H  is  the  total  layer  depth,  AT  = 

r  /  5  \  \  ~  r,-r2,  g  is  the  gravitational  acceleration,  a  is 

y  J  \  V  I  the  thermal  diffusivity,  and  u  is  the  kine- 

—  —  matic  viscosity.  Since  Ra,  is  always  defined 

105  A^x3  ■  5^-—  =g  as  positive  if  the  layer  is  unstably  stratified, 

whereas  Ra2  changes  sign,  Ra,  was  used  by 
-  -  Merker  et  al.  (1979)  to  describe  the  layer  sta- 

_  _  bility.  Figure  7  shows  the  critical  Rayleigh 

z  i  i  i  i  number  variation  with  T ,  in  the  case  of  con- 

-5  0  4  io  20  30  40  stant  wall  temperature  (i.e.  Tv  =  constant) 

with  the  density-temperature  relation  repre- 
Figure  8.  Critical  Rayleigh  number  Rac  sented  by  a  polynomial  of  order  p  =  5  and  n 

as  a  function  of  T ,  with  T2  as  a  param-  =  1  (number  of  Galerkin  terms).  Figure  8 

eter  for  qw  =  constant  (constant  wall  shows  the  case  of  constant  wall  heat  flux 

flux)  (after  Merker  et  at.  1979).  density  (i.e.  <?w  =  constant).  In  both  figures, 

it  can  be  noted  that  in  the  region  below  iso¬ 
therm  Ti  =  4°C,  the  water  layer  has  a  density  profile  with  no  maximum  value,  i.e.  the  layer 
is  completely  unstably  stratified.  However,  the  nonlinearity  of  the  density  profile  is  slight 
and  the  effect  on  the  critical  Rayleigh  number  is  not  large.  Rac  is  in  the  range  of  1708  to 
*  3600  for  Tv  =  constant  and  720to  *  1600  for  qv  =  constant.  On  the  other  hand,  for  the  re¬ 
gion  above  isotherm  r2  =  4°C,  the  layer  has  a  density  profile  containing  a  maximum  value 
between  the  two  boundaries  and  the  layer  is  only  partially  unstably  stratified.  The  nonlinear¬ 
ity  of  the  density  profile  is  strong  and  the  effect  on  Rac  is  great.  The  Rac  values  are  signifi¬ 
cantly  larger  than  those  obtained  from  the  classical  problem.  For  constant  wall  temperature 
and  rigid  boundary  conditions,  this  work  should  be  comparable  to  the  work  of  Sun  et  al. 


v 


(1969)  in  which  a  cubic  density-temperature  relation  was  used  and  the  entire  layer  depth  (as 
in  this  work)  was  used  in  defining  the  critical  Rayleigh  number. 

EXPERIMENTAL  STUDIES  OF  THE  ONSET  OF  CONVECTION  IN  A 
CIRCULAR  HORIZONTAL  MELT  LAYER 

The  first  experimental  study  that  focused  solely  on  determining  the  onset  of  convection  in 
a  water  layer  formed  by  melting  ice  was  reported  by  Yen  (1968).  In  a  study  of  the  effect  of 
buoyancy  on  the  melting  and  freezing  process,  Boger  and  Westwater  (1967)  claimed  the  criti¬ 
cal  Rayleigh  number,  Rac(Ra  =  gH}  AT0/av)  in  a  system  involving  phase  change  and  den¬ 
sity  inversion,  was  around  1700,  nearly  identical  to  the  value  reported  by  Globe  and  Dropkin 
(1959)  and  Schmidt  and  Silveston  (1959),  who  heated  the  horizontal  layers  of  water  from  be¬ 
low  but  did  not  test  in  the  region  of  the  density  inversion.  However,  to  have  Rac  *  1 700  as  in 
the  classical  case  of  normal  fluid  (i.e.  density  is  a  monotonic  function  of  temperature),  prop¬ 
er  values  of  H,  0,  A T,  a,  and  v  had  to  have  been  used  depending  on  the  direction  of  melting 
(upward  or  downward)  and  the  thermal  boundary  conditions.  If  the  ice  is  at  the  top  (i.e. 
melting  from  below),  AT  =  T- 4°C,  H  is  the  total  layer  depth,  0  is  evaluated  at  T„  and  the 
other  properties  of  a  and  v  are  evaluated  at  'A(T,  +  4°C).  If  the  ice  is  at  the  bottom  (i.e.  melt¬ 
ing  from  above),  they  proposed  using  AT  =  4°C,  H  is  the  portion  of  water  layer  depth  and 
was  evaluated  by  equating  the  heat  transported  through  the  convective  layer  to  the  heat  con¬ 
ducted  through  the  ice,  i.e.  H  =  -4 dxkt/k{T{,  in  which  d„  fcj,  and  7]  are  thickness,  thermal 
conductivity,  and  temperature  (maintained  at  the  end  of  the  sample  of  ice),  kt  is  the  effective 
thermal  conductivity  of  water  of  the  convective  layer,  0  is  evaluated  at  0°C,  and  a  and  v  are 
evaluated  at  the  mean  temperature  of  the  buoyant  region,  i.e.  2°C.  In  Yen’s  work,  a  rather 
large,  cylindrical,  bubble-free  ice  sample  was  used  (12.7  cm  in  diameter  and  25  cm  high)  and 
at  the  beginning  of  the  experiment  there  was  only  one  solid  phase.  A  constant  temperature 
boundary  was  applied  at  the  lower  end,  which  resulted  in  melting  from  below.  Because  of  the 
rapid  transition  of  the  melt  layer  from  the  conductive  to  the  convective  state,  it  restrained  the 
upper  limit  of  the  temperature  that  could  be  imposed.  For  T  ranging  from  7.72  to  25.2°C 
and  using  the  criteria  recommended  by  Boger  and  Westwater  (1967)  in  evaluating  Rac,  he 
found  on  the  contrary  that  the  Rac  values 
were  a  strong  function  of  T  and  can  be  rep-  'G«  io 5 
resented  by 

Rac  =  14,200 exp(-6. 64 x  10  !  7",).  (6) 

The  initial  temperature  of  ice  T0  was  found 
to  have  had  no  significant  effect  on  the  Rac 
values.  Figure  9  strongly  shows  the  variation 
of  Rac  with  the  lower  boundary  temperature. 

It  should  be  emphasized  that  the  study  of 
the  melt  layer  is  quite  distinct  from  the  clas¬ 
sical  problem  of  a  horizontal  layer  of  fluid 
of  invariant  depth.  For  a  water  layer  formed 
by  melting,  the  layer  depth  not  only  grows 
with  time  but  simultaneous  addition  of  water 
is  necessary  to  make  up  volume  shrinkage 
due  to  phase  transition  (an  intricate  device 
has  to  be  found  to  conduct  studies  of  this  na- 


a.  Melting  from  below. 


b.  Melting  from  above. 


Figure  10.  Schematic  representation  of  the  interde¬ 
pendence  of  stable  and  unstable  regions  with  the 
thermal  boundary  conditions. 

ture).  The  water  layer  can  be  formed  either  by  melting  from  below  or  above.  For  the  case  of 
melting  from  below  (i.e.  7,  >  4°C,  T2  =  0°C;  for  the  special  case  of  7,  =  4°C,  the  entire 
water  layer  will  always  be  stable),  an  unstably  stratified  region  will  extend  from  the  lower 
boundary  up  to  the  4°C  isotherm  and  a  stable  region  above  the  4°C  isotherm  will  extend  to 
the  upper  boundary  (the  water/ice  interface).  On  the  other  hand,  for  the  case  of  melting 
from  above  (for  7,  >  4°C,  7,  =  0°C),  there  will  be  a  stable  layer  adjacent  to  the  upper 
boundary  and  an  unstable  region  between  the  4°C  isotherm  and  the  water/ice  interface, 
while  for  T2  =  4°C,  the  entire  melt  layer  will  be  unstably  stratified.  Figure  10  shows  the  in¬ 
terdependence  of  the  stable  and  unstable  regions  with  the  thermal  boundary  conditions.  In  a 
follow-up  study  on  melting  from  below  and  above,  Yen  and  Galea  (1959)  reported  that  the 
critical  layer  depth  //c  (at  which  the  heat  transfer  mode  changes  from  conduction  to  convec¬ 
tion)  could  be  determined  either  by  locating  the  inflection  point  on  the  melting  front  vs  time 
or  by  locating  the  sudden  jump  of  the  temperature  gradient  in  the  stable  region.  Figure  1 1 
shows  the  variation  of  the  temperature  gradient  as  a  function  of  time  and  the  upper  boun¬ 
dary  temperature  (cases  of  melting  from  above).  It  can  be  seen  that  the  higher  the  upper 
boundary  temperature,  the  more  shallow  the  temperature  gradient  jump  at  the  onset  of  con¬ 
vection  Figure  12  compares  the  Rat  values  for  melting  from  below  and  from  above.  Experi¬ 
mental  Ra  values  were  evaluated  using  appropriate  values  of  critical  layer  depth  Hc:  A  = 

1 1  T  .  AT  r  for  melting  from  below  and  A  =  ^7  =  -T2  for  melting  from 

above  Theoretical  Ra  values  were  determined  from  a  graph  developed  by  Sun  et  al.  (1969) 
1 1  ig  *  >n  tho  paper  i  ,ii  er  evaluating  parameters  X,  and  \;  from  eq  4  and  5.  Figure  13  demon- 
oratc'  another  *av  ot  comparing  Ra,  values  using  7  or  7.  explicitly  as  a  variable.  This  fig- 
■  >e..>vered  I  ?tom  '  '2  to  2<  2  C  and  7  from  6.3  to  13. 09°C.  The  corresponding  ranges  of 
ii»:  an  ■  ' r om  I  'f’o  to  0  830  and  X; from  -0.31 1  to  0.390  for  melting  from  below, 

o.i  ■  >m  144410  l  M8  and  s  from  0  1 12  to  0.508  for  melting  from  above. 
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Figure  13.  Comparison  of  experimental 
and  theoretical  Rac  using  T,  or  T2  as  the 
variable  (after  Sun  et  al.  1969). 


Figure  14.  Comparison  of  Taf 
[2(1  -  /ifir4}  and  Rac/(2\,n •)  vs  X 
(after  Yen  1980). 


The  experimental  results  of  Yen  (1968)  and  Yen  and  Galea  (1969)  were  also  compared  with 
the  analytical  work  of  Veronis  (1963)  and  the  experimental  work  of  Legros  et  al.  (1974).  Ver- 
onis  considered  the  case  equivalent  to  melting  from  above  while  Legros  et  al.  studied  the  case 
equivalent  to  melting  from  below.  In  both  cases,  however,  the  layer  was  invariant,  while  in 
the  studies  of  Yen  (1968)  and  Yen  and  Galea  (1969)  the  water  layer  was  formed  (initially  at 
zero  depth)  continuously  as  a  resultant  of  phase  transition.  Veronis’  results  were  plotted  in 
terms  of  Tac/2(l-/i)'ir4  vs  Xin  which  Tac  is  the  critical  Taylor  number.  The  work  of  Legros  et 
al.  was  expressed  in  terms  of  Rac/2X‘ir4,  in  which  Ra  is  defined  as 

Ra  =  gy(AT)x  d\/ on 

where  y  is  a  temperature  coefficient  as  defined  in  eq  I,  AT  is  the  temperature  difference 
across  the  unstable  layer  and  in  the  cases  of  melting  from  above  and  below,  AT  =  4°C  and 
F, -4° C  respectively  [due  to  the  presence  of  (AT)2  in  the  definition  of  Ra,  this  comparison 
should  be  valid  for  both  cases  of  melting].  The  values  are  evaluated  as  follows:  For  melt¬ 
ing  from  above  dQ  =  -( 4/T2)Hc ,  and  for  melting  from  below  dc  =  [(7',-4)/7',J  Hc.  Figure  14 
shows  a  remarkable  agreement  among  these  studies, 

Some  unique  features  were  observed  in  the  studies  reported  by  Yen  (1968)  and  Yen  and 
Galea  (1969)  during  the  formation  and  continuous  deepening  of  the  water  layer.  Prior  to  the 
onset  of  convection,  the  water/ice  interface  remains  a  planar  surface  and  the  prevailing  heat 
transfer  mode  is  conduction.  Cells  start  to  form  as  soon  as  the  Rac  value  is  attained.  At  the 
very  beginning  stage,  these  cells  possess  regular  patterns  and  are  hemispherical  and  higher  in 
the  center.  They  have  a  circular  cross-section,  differing  from  the  well-known  hexagonal  B£- 
nard  cells,  and  are  distributed  evenly  over  the  entire  water/ice  interface.  Figure  15a  shows 
schematically  the  interface  morphology  observed  shortly  after  the  onset  of  convection  in  the 
case  of  melting  from  below.  It  is  interesting  to  note  that  in  the  case  of  melting  from  above,  a 
stable  region  would  form  over  the  unstable  region  adjacent  to  the  interface  as  soon  as  the 
melted  layer  was  formed.  The  interface  (see  Fig,  15b),  in  contrast  to  the  case  of  melting  from 
below,  is  covered  with  a  series  of  circular  concentric  ridges  equally  spaced  from  each  other. 
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a.  Melting  from  below.  b.  Melting  from  above. 

Figure  15.  Schematic  of  water/ice  interface  shortly  after  onset  of  convection  (after  Yen  and  Galea 
1969). 


The  height  between  crest  and  trough  increases  with  time  and  it  can  reach  a  value  as  large  as 
5.0  mm.  However,  this  pattern  was  replaced  by  that  of  small  inverted  hemispherical  cells  that 
gradually  enlarged  and  finally  formed  an  irregular  and  unsymmetrical  interface  as  the  melt¬ 
ing  progressed  (as  in  the  case  of  melting  from  below). 


TEMPERATURE  STRUCTURE  AND  HEAT  TRANSFER 


In  a  horizontal  layer 

Townsend  (1964)  was  the  first  to  make  a  detailed  study  of  the  temperature  structure  and 
natural  convective  heat  transfer  in  water  over  an  ice  surface.  The  experimental  tank  consist¬ 
ed  of  a  30-  x  30-cm  bottom  made  of  Dural  0.95  cm  thick  with  sides  made  of  Perspex  of  the 
same  thickness.  Distilled  water  filled  the  tank  to  a  depth  of  15  cm,  and  the  tank  bottom  was 
cooled  by  a  shallow  dish  filled  with  liquid  nitrogen  (not  in  direct  contact  with  the  bottom, 
however;  it  should  be  noted  that  the  purpose  of  bottom  cooling  was  only  to  maintain  the 
lower  boundary  at  0°C).  The  free  water  surface  temperature  was  maintained  electrically 
through  another  Dural  plate  suspended  about  5  mm  above.  The  most  striking  feature  of  the 
steady  temperature  distribution  was  the  creation  of  a  region  of  constant  mean  temperature. 
Townsend  indicated  that,  after  establishing  a  constant  mean  temperature  region,  the  only 
observed  change  was  the  upward  extension  of  its  boundary.  Figure  16  shows  the  temperature 
distributions  during  the  approach  to  equilibrium.  Figure  17a  shows  the  equilibrium  distribu¬ 
tion  of  the  mean  temperature.  It  can  be  seen  that  the  constant  temperature  region  has  a  mean 
temperature  of  3.2°C,  significantly  below  the  temperature  of  maximum  density  (i.e.  =  4°C). 
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Figure  16.  Temperature  distribution  during 
the  approach  to  equilibrium.  The  numbers  by 
the  curves  are  the  elapsed  times  in  minutes  from  the 
start  of  cooling  (after  Townsend  1964). 
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b.  Stable  and  disturbed  regions 
(note  scatter  of  points  for  tempera¬ 
tures  near  6  °C,  indicating  unusual¬ 
ly  large  fluctuations). 


Figure  1 7.  Mean  temperature  distribution  for  thermal  equilibrium  corrected  to  a 
surface  temperature  of  23  °C  (after  Townsend  1964). 


Townsend  also  noted  the  peculiarities  of  “overshoot”  regions  at  each  end  where  the  gradient 
of  the  mean  temperature  changes  sign  (Fig.  17a).  The  mean  temperature  increases  with 
height  and  reaches  3.6°C  before  decreasing  to  3.2°C;  at  the  upper  limit  of  the  region  it  falls 
to  2.9 °C  before  increasing  rapidly  in  the  stable  region  above.  The  upper  overshoot  was  post¬ 
ulated  to  be  a  consequence  of  the  slowing  down  of  the  lower  parts  of  rising  columns  of  cold 
water  as  their  tops  enter  the  stable  layer.  However,  Townsend  stated  that  the  lower  over¬ 
shoot  cannot  be  interpreted  in  a  similar  manner,  as  descending  columns  of  hotter  water  are 
not  observed,  and  he  indicated  this  may  be  attributed  to  the  subsidence  of  water  with  a  tem¬ 
perature  between  3.2  to  4°C.  Townsend  also  noted  large  fluctuations  of  temperature  just 
within  the  region  of  stable  stratification  but  they  fell  off  very  rapidly  with  height  and  were 
very  small  if  the  local  temperature  exceeded  10°C  (Fig.  17b). 

Townsend  evaluated  the  downward  heat  flux  q  (i.e.  the  heat  extracted  through  the  tank 
bottom)  by  fixing  z  =  zr  as  a  reference  plane  in  the  stable  region  and  assuming  horizontal  ho¬ 
mogeneity  as  follows: 


1 1 


(7) 


-  i.  — »  41.  4V*'A'«  tl-  a|  ..I  k.l  <(i  ».t  >-■  » ,1 


I 


c*www 


■-^4  i  <r-rr>*  +  *(||)i 

A  \  /  I 


where  the  first  term  on  the  right  represents  the  rate  of  change  of  heat  content  below  fixed 
plane  z„  T is  the  local  mean  temperature,  Tt  is  the  reference  temperature  (taken  to  be  3.2  °C), 
q,  cp,  and  k  are  the  density,  heat  capacity,  and  thermal  conductivity  of  the  water,  and  r  is  the 
distance  measured  vertically  upward  from  the  ice  surface. 

With  the  measurement  of  T  both  in  the  stable  and  unstable  regions  as  a  function  of  time, 
the  two  terms  on  the  right  of  eq  7  can  be  evaluated.  The  heat  flux  was  reported  to  be  from  34 
mW/cnrJ  in  the  beginning,  decreasing  to  about  26  mW  cm‘J  as  the  temperature  distribution 
approaches  equilibrium.  He  attributed  this  reduction  to  heat  intake  through  the  side  walls  of 
the  test  tank  as  the  thick,  cold,  constant-temperature  region  formed.  However,  since  there 
were  no  changes  in  temperature  distribution  and  the  pattern  of  convection,  he  concluded 
that  the  equilibrium  heat  flux  in  a  horizontally  homogeneous  system  would  be  close  to  34 
mW  cm  '.  Using  this  flux,  Townsend  derived  the  following  expression: 

Q  =  2-  =  0. 1 56(AT)i/}  (yiig)w)  (a/v)in  (8) 

QCp 

in  which  -47~is  the  temperature  difference  between  the  tank  bottom  and  the  temperature  of 
maximum  density  Tmax  =  3.98°C  and  y  is  defined  in  eq  1.  Equation  8  is  the  dimensionally 
analogous  form  for  the  heat  loss  from  a  smooth  heated  plane  involving  fluids  with  linear  ex¬ 
pansion  given  by 

Q  =  0.193(A7)4/3(a/jg)1/3(tt/e)|/3  (9) 

in  which  AT  is  the  temperature  difference  between  the  surface  and  far  above  the  plane. 

Myrup  et  al.  (1970)  also  studied  this  convection  problem  over  ice,  which  they  found  has 
great  value  in  understanding  and  modeling  geophysical  and  astrophysical  convection  sys¬ 
tems.  Their  experimental  set-up  was  similar  to  Townsend’s.  The  water  was  initially  at  room 
temperature  and  the  tank  bottom  was  maintained  slightly  higher  than  0°C  by  a  water  and  ice 
bath  beneath  it.  The  top  of  the  tank  was  covered  with  an  aluminum  plate  in  direct  contact 
with  the  water  surface  that  was  conditioned  to  equilibrium  with  the  room  air.  They  made  ex¬ 
tensive  measurements  on  mean  temperature  profiles  as  well  as  temperature  fluctuations  at 
every  level.  Temperature  profiles  similar  to  those  of  Townsend  were  reported.  Figure  18 
shows  a  reproduction  of  the  temperature  record  at  a  single  height  of  2.8  cm  above  the  lower 
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Figure  18.  Temperature  at  2.8  cm  above  the  tower  boundary  maintained  near 
O'C  (after  Myrup  et  at.  1970). 
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plate  (bottom).  It  is  evident  that  distinctly  different  kinds  of  temperature  fluctuations  oc¬ 
curred.  At  about  160  to  220  minutes  after  the  beginning  of  the  experiment  there  were  some 
sort  of  periodic  warm  fluctuations,  some  of  which  were  warmer  than  4°C.  This  phenomenon 
was  attributed  to  the  effects  of  gravity  waves  bobbing  up  and  down  in  the  stable  layer  just 
above  the  convective  layer.  From  220  minutes  on,  the  temperature  fluctuations  distinctly 
changed.  The  irregular  cold  deviations  are  due  to  penetrations  by  convective  plumes  of  the 
buoyant  cold  fluid  rising  from  the  lower  boundary  layer.  Photographs  delineated  the  forma¬ 
tion  and  subsequent  development  of  three  distinct  layers.  The  boundary  layer,  nearest  the 
lower  plate,  is  occupied  by  the  densest  concentration  of  the  dye;  the  dark  upper  portion  is  the 
fluid  warmer  than  4°C  into  which  the  dye-laden  convective  currents  have  not  penetrated, 
and  in  between  is  the  convective  layer.  Myrup  et  al.  (1970)  made  no  attempt  in  their  study  to 
evaluate  the  heat  fluxes. 

In  a  circular  horizontal  melt  layer 

Along  with  the  study  of  the  determination  of  the  onset  of  convection  in  the  melt  layer,  the 
mean  temperature  of  the  growing  melt  layer  was  measured  intermittently  (Yen  1 968,  Yen  and 
Galea  1969).  The  ice/water  system  in  these  studies  was  drastically  different  from  those  in  the 
works  reported  by  Townsend  (1964)  and  Myrup  et  al.  (1970)  in  that  only  solid  ice  was  present 
at  the  onset  of  the  experiment.  The  depth  of  the  melt  layer  started  from  zero  and  deepened 
continuously  as  the  melting  progressed  either  from  above  or  below.  Thermistors  with  a 
0.02-s  response  time  were  used  to  measure  intermittently  the  temperature  profile  of  the  grow¬ 
ing  layer.  To  minimize  possible  disturbance  of  the  temperature  field  when  the  thermistors 
were  inserted,  they  were  installed  at  different  radii  and  distributed  over  the  circular  area. 
Typical  mean  temperature  profiles  are  shown  in  Figure  19,  for  melting  from  below,  and  in 
the  following  figure,  for  melting  from  above.  Each  data  point  represents  four  thermistor 
readings  from  the  same  layer  depth.  It  is  noted  that  for  both  melting  cases  the  thermistor 
readings  contained  a  random  component  associated  with  the  fluctuating  temperature  fields 
in  the  convection  region,  demonstrating  that  the  convective  motion  was  unsteady.  Figure  19a 
shows  the  temperature  profile  just  before  the  onset  of  convection  for  71  at  4.02,  6.65,  7.46, 
8.11,  and  9.84°C  respectively,  with  initial  temperature  roat  -11.5°C.  It  dearly  shows  that, 
for  4.02 °C  at  418  minutes  after  the  beginning  of  the  experiment,  the  temperature  profile  re¬ 
mains  linear  as  expected  (theoretically  the  melt  layer  will  remain  stable  indefinitely  as  long  as 
71  »  4°C).  However,  for  71  at  9.84°C,  the  convection  was  about  to  commence  barely  36 
minutes  after  the  experiment  started.  This  demonstrated  strongly  the  dependence  of  the  layer 
stability  on  the  imposed  71 .  At  the  onset  of  convection  (Fig.  19b),  the  temperature  curves  are 
all  deviated  upward  from  the  linear.  In  all  these  figures,  the  number  by  each  curve  is  the  time 
elapsed  from  the  initiation  of  the  experiment;  the  times  chosen  were  nearly  equal  to  each 
other  to  show  the  significance  of  71  affecting  the  onset,  and  implicitly  the  intensity,  of  the 
convective  mixing  of  the  unstable  region.  Figure  19c  shows  the  pseudo-steady-state  tempera¬ 
ture  distribution  of  the  melt  layer;  it  consisted  of  a  nearly  constant  temperature  zone,  but 
this  constant  temperature  was  noted  to  be  dependent  on  the  71  imposed,  i.e.  the  higher  the 
value  of  71,  the  higher  the  temperature  of  the  constant  temperature  zone.  The  four  thermis¬ 
tor  readings  at  each  level  exhibited  considerable  fluctuations,  especially  at  the  intersection 
between  the  convective  and  the  upper  conductive  layer  adjacent  to  the  water/ice  interface. 

Figure  20  shows  some  typical  mean  temperature  profiles  of  the  melt  layer  formed  by  melt¬ 
ing  from  above.  Contrary  to  the  previous  cases,  the  layer  is  most  stable  when  71  is  the  highest 
(because  the  driving  force  to  create  the  convective  motion  merely  results  from  the  density  dif¬ 
ferences  between  0°  and  =  4°C).  Figure  20a  shows  the  mean  temperature  profiles  just  before 
convective  motion  starts  for  7j  ranging  from  1I.92°C  to  as  high  as  39.75 "C '.  It  can  be  ob¬ 
served  that  the  melt  layer  became  progressively  more  stable  as  71  increased,  f  igure  20b  dis- 
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b.  Development  of  the  convective  layer. 


c.  Pseudosteady  state. 

Figure  19.  Mean  temperature  profiles  for 
melting  from  below  (after  Yen  1984). 


tinctly  shows  that  (after  almost  the  same  time  span),  for  T7  =  39.75 °C  and  at  t  =  238  min, 
the  temperature  is  just  beginning  to  deviate  from  the  linear  distribution,  but  on  the  other 
hand,  for  Tz  =  11.95°C  and  at  t  =  256  min,  a  constant  temperature  region  is  already  well 
developed.  Figure  20c  shows  clearly  the  establishment  of  the  well-developed  constant  tem¬ 
perature  zone  as  shown  in  Figure  19c  but,  contrary  to  the  case  of  melting  from  below,  the 
constant  temperature  zone  is  =  3.2 °C  regardless  of  what  value  of  T:  was  imposed  on  the  up¬ 
per  boundary.  Slight  temperature  reversals  can  be  observed  in  Figures  19c  and  20c.  Over¬ 
shoots  were  noted  in  all  temperature  traverses  (regardless  of  whether  the  melting  was  from 
above  or  below)  in  the  pseudosteady  conditions,  substantiating  the  findings  of  Townsend 
(1964). 
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Figure  21.  Nondimensional  mean  temperature  profile  for  melt¬ 
ing  from  above  (after  Yen  1980). 

The  molecular  scales  are  based  on  the  assumption  that  in  certain  regions  of  the  flow  the 
convection-layer  depth  is  not  dynamically  significant,  but  molecular  diffusion  is  important. 
Hence,  the  molecular  scales  are  expected  to  be  appropriate  in  the  boundary  layer  regions 
such  as  the  conduction  layer  and  the  interfacial  layer.  On  the  other  hand  the  convection 
scales  are  intended  to  apply  in  the  core  of  the  convection  layer  where  the  length  scale  of  the 
convection  motion  is  on  the  order  of  the  layer  depth.  Figure  21  shows  the  typical  nondimen¬ 
sional  mean  temperature  profile  for  the  case  of  melting  from  above,  including  the  works  of 
Townsend  (1964)  and  Adrian  (1975)  on  convection  of  water  over  ice  (a  thin  layer  at  maxi¬ 
mum).  It  should  be  pointed  out  that  in  Townsend’s  and  Adrian’s  work,  the  water  layer  did 
not  result  from  melting;  instead  an  invariant  and  rather  deep  water  layer  (=15  cm)  was  used. 
Adrian’s  results  seem  to  agree  remarkably  well  with  Yen’s  (1984). 

Figure  22  shows  the  nondimensional  mean  temperature  profile  for  melting  from  below. 
The  shape  of  the  curve  is  the  opposite  of  the  curve  shown  in  Figure  21 .  It  should  be  pointed 
out  that  the  temperature  of  the  convection  layer,  indicated  in  Figure  !9c  as  depending  on  the 
value  of  r,,  the  nondimensional  mean  temperature— derived  from  the  limited  available 
data— fell  more  or  less  on  a  single  curve  as  shown  in  Figure  21 ,  with  a  somewhat  higher  ratio 
of  H/0o  for  the  convection  layer.  The  only  experimental  work  closely  related  to  the  case  of 
melting  from  below  was  conducted  by  Legros  et  al.  (1974)  for  those  experiments  with  an  up¬ 
per  boundary  less  than  4°C.  However,  their  work  was  aimed  at  determining  the  critical  tem¬ 
perature  difference  across  the  layer,  and  no  attempt  was  made  to  measure  the  temperature 
profile  within  the  layer. 

Yen  (1980)  also  derived  heat  flux  expressions  for  the  melt  layer  formed  by  melting  from 
both  below  and  above.  For  experiments  in  which  the  melting  rates  were  determined,  the  total 
upward  or  downward  heat  flux  for  melting  from  below  or  above  was  evaluated  by 
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The  term  (QCp/2t)[z(t)(T +  Tm)}  represents  the 
mean  sensible  heat  content  variation  of  the  en¬ 
tire  layer  of  a  depth  z(t).  In  the  case  of  melting 
ice,  Tm  =  0°C  and  T  stands  for  either  T ,  or  7). 
The  contribution  of  this  term  to  the  overall  heat 
flux  was  found  to  be  much  more  significant  in 
the  case  of  melting  from  above.  For  those  ex¬ 
periments  in  which  the  mean  layer  temperature 
was  measured  as  the  melting  progressed,  the 
heat  flux  was  approximated  by 


\(d_T\  (dT)\, 

l  '  dzn  +  dzn  I 


u)  9,6  -V.  in  which  (dT/dz)i  and  (dT/dz)\  are  the  mean 

i _  .  _ ;  temperature  gradients  at  the  stable  region  near 

0  5  1 0 

g/g  the  upper  boundary  for  melting  from  above 

and  the  lower  warm  plate  for  melting  from  be- 
Figure22.  Nondimensional mean  tern-  ,ow  Subscripts  ,  and  2  indicate  the  beginning 

perature  profile  for  melting  from  be-  ,  ,  .  ,  .......  .  . 

<  ,  and  end  of  each  period  and  At  is  the  total  time 

low  I after  Yen  1980).  .  ,  „  ,  , 

period.  For  calculating  a  value  of  q,  at  least  a 

dozen  or  more  of  these  periods  (with  varying 

durations)  were  used.  The  heat  flux  from  above 

was  found  to  be  a  weak  function  of  7)  and  can  be  represented  by 


q  =  177(7))° 303 


with  an  average  q  -  474  W/mJ.  The  above  expression  is  valid  for  7)  ranging  from  1 1.75  to 
39.9°C.  For  melting  from  below.  Yen  (1980)  reported 


q  =  -1900  +  3 1 5  ( 7~, ).  (15) 

For  T,  ranging  from  7.7  to  25.5°C,  the  heat  flux  is  found  to  be  a  linear  function  of  7) .  High¬ 
er  values  of  T,t  and  thus  higher  temperatures  in  the  convective  zone,  resulted  in  greater  con¬ 
vective  motion  in  the  unstable  region  and  consequently  reduced  the  thickness  of  the  stable 
layer  adjacent  to  both  the  ice  and  the  thermal  boundary.  The  works  of  Townsend  (1964)  and 
Adrian  (1975)  are  similar  to  the  case  of  melting  from  above.  However,  in  their  investiga¬ 
tions,  a  rather  deep  invariant  water  layer  was  used  throughout  the  experiment,  and  there  was 
no  phase  transition  taking  place  at  the  bottom  of  the  tank.  They  both  reported  a  nearly  iden¬ 
tical  heat  flux  of  approximately  340  W/m!  independent  of  the  initial  water  temperature.  The 
discrepancy  in  reported  q  values  (in  Yen’s  work  the  melt  layer  was  formed  as  melting  pro¬ 
gressed,  and  in  Townsend’s  and  Adrian’s  case  there  was  only  a  liquid  phase  to  begin  with) 
are  believed  to  be  due  to  the  real  difference  in  the  transfer  processes  involved  in  the  experi¬ 
mental  system. 


HEAT  TRANSFER  STUDIES  IN  NONPLANAR  GEOMETRIES 


The  most  comprehensive  and  earliest  theoretical  study  of  melting  free  convective  heat 
transfer,  other  than  the  horizontal  geometries  such  as  spheres  and  cylinders,  was  reported  by 


Merk  (1954).  Employing  a  third-order  density-temperature  polynomial  of  water  density  and 
applying  the  Von  Karman-Pohlhausen  integral  method  for  the  case  of  Pr  >>  1,  he  success¬ 
fully  solved  the  boundary  layer  equation  and  developed  a  general  Nusselt  number  ratio: 

6715  _  8471  „  \[1  +(S/2)]311/4 


where  S  is  the  shape  factor  of  the  temperature  profile  and  is  connected  to  the  Stefan  number 
St  =  cp(Tm-TJ/L  as 

s~-2*Wt  -  I  V'-TSl 


^  "  (S  +  2)J  ’ 

where  Tm  is  the  bulk  water  temperature.  The  values  of  P,,  Pi,  and  P,  are  expressed  as 

p  -  i  _  ^  5  +  _L2_ 

P'  ~  1  217  5  +  434  A 


p  =  ,  i21  c:_  J9_5, 

1343  2686  2686 


P  _  6475  3537  199  19 

’  8471  16942  8471  ^  +  16942 


The  values  of  m  and  rt  are  defined  respectively  as 


m  =  (0i  +  30,T  )6m/N@ 


n  =  M'JWPJ 


where  &m  =  Tm-  T  and  N  and  0  are  defined  as 


N=  \+0lTm  +  l3iT>  +  0,T> 


0„=W,+2l5iT„  +  3l3,T3/N  (25 

in  which  the  0s  are  the  coefficient  in  the  formula  for  the  specific  volume  of  the  water,  i.e. 


—  =  —  (1  +  0,r+  0iT‘  +0,T'). 

e  Co 


Nun  is  the  Nusselt  number  for  the  case  of  St  =  m  =  n  =  0.  Merk  (1954)  reported  that  for 
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Figure  23.  The  influence  of  melting 
(S>  <  0)  and  solidification  (S  >  0) 
on  the  heat  transfer  in  thermal  con¬ 
vection  for  large  values  of  the 
Prandtl  number  (after  Merk  1954). 


large  values  of  the  Prandtl  number,  neither  the  shape  of  the  body  nor  the  position  on  the  sur¬ 
face  influences  the  ratio  Nu/Nu0.  Therefore  this  ratio  can  be  replaced  by  Nu/Nu0,  and  can 
be  expressed  by  the  following  if  only  the  effect  of  melting  is  considered  (i.e.  for  m  =  n  =  0) 


Figure  23  is  a  graphic  representation  of  eq  27;  it  clearly  shows  that  for  S  <  0  or  for  melting, 
Nu  <  Nu0,  while  for  solidification,  i.e.  S  >  0,  Nu  >  Nu0.  For  the  case  of  no  melting  but 
with  convective  inversion,  then  S  =  0  and  hence  =  P2  =  P,  =  1 .  Equation  26  becomes 

Nu 

~  =  (1  +  0.4912m  +  0.2730n)l/4.  (28) 

Nu0 

By  defining  Nu0  for  large  values  of  the  Prandtl  number  as 

Nu0  =  C(GrPr)1/4  =  0^".  (29) 
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Equation  28  can  be  written  as 


10  Jm(l  +  0.491 2m  +  0.2730n)J1/4 


where  /  is  a  dimensionless  number.  For  small  values  of  T  ,  0 

oo  T  ^  oo 

0X,  and  n  =  0,  eq  30  can  be  simplified  to 


(30) 


0 ,  +  2 0tT.  m  *  0>9m/ 


*,*  V  v « 
«  #  *  * 

*  “w  *,  I 


/=  (1. 509 A  Tm-T„  Tm-Tm  ),/4 


(31) 


8  12 
T„  it) 


Figure  24.  Convective  inversion  for  melt¬ 
ing  ice  in  water  of  temperature  T  a>.  The  ar¬ 
rows  indicate  the  direction  of  the  flow  along 
the  surface.  The  dashed-dotted  curve  shows 
the  behavior  neglecting  melting,  while  in  the 
full  curve  the  melting  is  taken  into  account 
(after  Merk  1954). 


in  which  Tw  is  the  inversion  temperature: 

7]v  =  -0.663  §•  -  0.326 Tm.  (32) 

P2 

Using  appropriate  values  of  0,  and  02  and  taking  Tm  =  0°C,  Merk  reported  a  value  of  7jv  = 
5.005 °C,  indicating  that  the  inversion  temperature  for  melting  ice  in  water  is  somewhat 
greater  than  =  4°C.  The  significance  of  the  inversion  temperature  is  clearly  seen  from  eq  31 
at  Tn  =  Tiv,  Nu  =  0,  and  the  direction  of  the  flow  in  the  boundary  layer  along  the  surface  of 
the  body  is  inverted  (for  <  Tn  the  flow  is  upward  and  for  T x  >  7jv  it  is  downward).  Us¬ 
ing  the  general  eq  16,  Merk  derived  the  minimum  Nusselt  number  at  7]v  =  5.30°C  for  S  =  0 
(no  melting)  and  7jv  =  5.3  PC  for  S  <  0  (with  melting).  Figure  24  shows  the  effect  on  the 
value  of/(i.e.  the  value  of  Nu)  with  and  without  melting  with  the  convective  inversion.  It  is 
clearly  shown  that  the  effect  of  melting  is  only  appreciable  for  Tm  >  Tw  and  may  be  neglect¬ 
ed  for  Tx  <  7jv.  Experimental  results  of  Dumori  et  al.  (1953)  and  analytical  results  from  the 
nonmelting  vertical  plate  study  by  Ede  (1955)  generally  confirmed  Merk’s  findings. 

Tkachev  (1953),  using  photographic  techniques,  reported  a  minimum  Nusselt  number  for 
melting  ice  cylinders  at  5.5 °C  and  was  the  first  to  notice  the  peculiar  nature  of  the  maximum 
density  boundary  layer.  He  suggested  that  under  certain  conditions  the  boundary  layer  might 
be  split  with  a  predominantly  upward  motion  immediately  adjacent  to  the  ice  surface  and  a 
region  of  downward  motion  outside  this.  Tkachev  conducted  melting  experiments  on  spheres 
as  well  as  on  vertical  and  horizontal  cylinders.  Using  the  same  initial  cylinder  diameters  but 
with  various  bulk  water  temperatures  he  found  that  the  coefficient  of  heat  transfer  is  lowest 
for  a  water  temperature  of  about  5.5 °C.  He  correlated  his  data  with  the  following  dimen¬ 
sionless  expressions  as 

Numd  =  0.40(GrPrC  (33) 


Numd  =  0.104(GrPr)^d  (34) 

for  cylinders  for  the  values  in  the  range  of  10!  <  GrPr  <  107  and  (GrPr)md  >  107  respective¬ 
ly.  The  corresponding  equations  for  spheres  are 

Numd  =  0.54(GrPr)md  (35) 
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for  laminar  flow  [103  <  (GrPr)md  <  107]  and 

Numd  =  0.135(GrPr)"d4  (36) 

for  turbulent  motion  [(GrPr)md  >  10’].  Subscript  md  represents  the  physical  properties,  and 
the  diameter  of  the  sphere  or  cylinder  was  evaluated  at  the  arithmetic  mean  temperature  (i.e. 
( Tm  +  TX)/2J  and  diameter  [i.e.  (d0  +  df)/2]  where  d0  and  df  are  the  initial  and  final  diameter 
at  the  end  of  the  experiment.  Based  on  his  experimental  data,  Tkachev  further  presented  an 
expression  for  determining  the  time  required  for  completing  melting  as 

StFomdNumd  =  0.305.  (37) 


The  suggestion  of  split  boundary  layer  flow  was  verified  by  the  analytical  and  experiment¬ 
al  work  of  Schecter  and  Isbin  (1958),  in  which  an  isothermal,  vertical,  nonmelting  plate  was 
used.  Figures  25  and  26  compare  the  theoretical  and  experimental  results  for  the  unidirec- 


tional  and  inverted  convections  respectively.  The  theoretical  curves  in  Figures  25  and  26  are 
given  by 


Nu  =  0.892 


Gr 


/  1  m  n-J 

(t  +  T  +  T 


(0.952  +  Pr) 


Pr 


1'2 


(38) 


and 


Nu  =  0.652 


GrPr(|  +  f  * 


c 


(39) 


where  m  and  n  ate  defined  in  eq  22  and  23.  They  reported  that  a  test  of  the  type  of  region 
that  will  prevail  for  given  conditions  of  plate  and  bulk  temperature  can  be  stated,  i.e.  if 

1  m  n  _ 

y  +  y  +  y  S  0, 

there  will  be  normal  convection  (unidirectional),  and  if 


1  m  n 

T  +  T  +  ~  <  °’ 


there  will  be  inverted  convection.  They  concluded  that  the  heat  transfer  coefficient  can  be 
predicted  for  both  regions  with  a  deviation  in  Nusselt  number  of  ±  lO^o  provided  the  abso¬ 
lute  value  of  (1/3  +  m/5  +  n/7)  >  0.05  by  use  of  eq  38  or  39,  depending  on  the  convective  re¬ 
gion.  However,  it  should  be  noted  that  the  boundary  layer  equations  as  approached  either  by 
the  Von  Karman-Pohlhausen  integral  method  or  by  the  similarity  transformation  method 
did  not  yield  meaningful  results  under  split-flow  conditions. 

To  resolve  this  problem  Vanier  and  Tien  (1968)  used  an  accurate  numerical  method  to 
solve  the  similarity  equations  for  a  semi-infinite  vertical  plate  at  constant  temperature  Tw  im¬ 
mersed  in  an  indefinitely  large  volume  of  water  at  bulk  temperature  Tx.  They  reported  a  new 
solution  is  necessary  for  every  combination  of  7"w  and  Tx.  By  obtaining  several  hundred  such 
solutions,  the  authors  were  able  to  map  out  temperature  zones  for  each  flow  regime.  The 
split  boundary  layer  was  found  to  be  confined  to  two  distinct  triangular  regions  within  which 
the  similarity  equations  become  quite  intractable.  They  confirmed  the  findings  of  Merk 
(1954)  that  the  melting  heat  transfer  rates  were  closely  similar  to  the  case  of  nonmelting. 

Vanier  and  Tien  (1970)  conducted  experimental  work  aimed  at  relating  their  numerical 
plate  results  to  a  more  practical  geometry  of  the  sphere  (including  the  effect  of  changing 
body  configuration).  This  was  also  partially  motivated  by  lack  of  detailed  analysis  and  corre¬ 
lations  of  the  experimental  results  on  the  melting  of  ice  spheres  and  cylinders  presented  by 
Dumor£  et  al.  (1953)  and  Tkachev  (1953).  They  presented  their  results  in  a  least-squares- 
fitted  semiempirical  equation  as 

Nu  =  2  +  C(GrPr)14  (40) 

and  found  that  for  T m  >  7°C,  it  appears  that  the  results  are  not  affected  by  the  maximum 
density  and  the  best  value  of  C  is  0.422  ±  0.006  in  the  range  of  1 .7  x  10‘  <  GrPr  <  2.4  x  10* 
(Fig.  27,  curve  b).  However,  for  Tx  <  7°C,  nearly  the  same  value  of  Cis  found  but  with  con¬ 
siderably  more  scattering  (twice  the  standard  deviation),  indicating  the  need  at  least  to  incor- 
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Figure  27.  Susselt  number  as  function  of  Rayleigh  number 
(GrPr )  t after  l  amer  and  Tien  1970). 


Figure  28.  Susselt  number  variation  with 
bulk  ’emperature  for  approximately  con¬ 
stant  sphere  diameters  t after  l  amer  and 
Tien  1970). 


porate  another  parameter  such  as  T ^  to  describe  adequately  the  heat  transfer  under  these 
conditions.  This  was  done  by  separating  the  low  temperature  data  into  positive  and  negative 
deviations.  For  Tx  <  3.8°C,  Nu  values  are  higher  than  expected  (curve  a),  while  for  4.1 0  < 
T ^  <7.1  °C,  Nu  values  are  too  low  (curve  c).  To  check  the  one-quarter  power  assumption  in 
eq  40,  a  two-parameter  fit  was  carried  out,  resulting  in 

Nu  =  2  +  0.437(GrPr)0  24\  (41) 

which  provides  an  excellent  verification.  However,  if  theGrashof  number  were  calculated  by 
using  an  arithmetic  mean  temperature  basis  of  Tx,  the  constant  C  was  found  to  be  0.52  for 
Tx  >  14°C.  This  provides  a  remarkable  agreement  with  the  results  reported  by  Tkachev 
(1953).  The  effect  of  sphere  diameter  and  maximum  density  on  heat  transfer  can  be  seen  in 
Figure  28.  These  curves  are  in  general  agreement  with  the  flat  plate  results  reported  by  Van- 
ier  and  Tien  (1968),  which  show  a  sharp  minimum  between  5°  <  Tx  <  6°C.  To  ascertain  the 
effect  of  sphere  diameter,  Vanier  and  Tien  (1970)  proposed  a  correlation  of  the  sphere  results 
with  those  from  theoretical  analysis  of  a  melting  plate  by 

Ntip/Nu  =  C(L/D)1  *  (42) 

where  L  and  D  are  the  characteristic  height  of  the  plate  and  the  diameter  of  the  sphere.  The 
least-squares-fitted  constant  C  was  found  to  be  1.106  t  0.144.  The  scaled-up  experimental 
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data  are  shown  in  Figure  29.  This  figure  clearly  means  that  a  melting  sphere  behaves  very 
similarly  to  a  melting  flat  plate  and  that  if  all  the  transfer  parameters  are  equal  (including 
temperature,  characteristic  length,  and  surface  area),  about  11%  more  heat  is  transferred  to 
the  plate  than  to  the  sphere.  This  is  the  effect  of  curvature  on  the  flow  velocities  and  is  in 
good  agreement  with  the  analytical  results  reported  by  Merk  and  Prins  (1954)  for  nonmelting 
free  convection  systems  without  maximum  density  effects,  i.e. 


Nup/Nu  =  1.14  (L/D)3'4.  (43) 

The  minimum  Nusselt  number  for  spheres  occurs  at  Tx  =  5.35  ±  0.2 °C,  as  compared  to  the 
value  of  5.31  °C  based  on  Merk’s  (1954)  theoretical  results. 

The  most  recent  study  of  heat  transfer  and  ice  melting  in  ambient  water  near  its  density  ex¬ 
tremum  was  reported  by  Bendell  and  Gebhart  (1976).  In  their  experiment,  a  vertical  ice  slab 
(30.3  cm  high,  14.8  cm  wide,  3  cm  thick  initially)  was  immersed  in  water  at  a  uniform  bulk 
ambient  temperature,  T w.  Figure  30  shows  the  experimental  results  along  with  the  analytical 
results  of  Gebhart  and  Mollendorf  (1978)  and  the  results  predicted  with  the  Boussinesq  ap¬ 
proximation.  Gebhart  and  Mollendorf  s  work  is  similar  to  that  of  Vanier  and  Tien  (1968)  ex¬ 
cept  it  gives  a  more  accurate  representation  of  the  density-temperature  of  water.  As  Vanier 
and  Tien  pointed  out,  in  the  inversion  region  the  validity  of  the  simplest  boundary  layer 
theory  becomes  questionable.  However,  beyond  that  region,  the  experimental  Nusselt  num¬ 
ber  values  were  nearly  equal  to  those  predicted  by  theoretical  analysis.  Bendell  and  Gebhart 
(1976)  reported  that  for  a  melting  vertical  ice  surface,  upflow  occurred  when  Tx  s  5.6°C. 
For  T x  >  5.5°C,  downflow  was  observed  and  was  found  to  be  in  good  agreement  with 
earlier  results.  They  found  the  minimum  Nusselt  number  for  the  experimental  temperature 
range  2.2°  <  Tx  <  25.2°C  occurred  at  Tx  =  5.6°C.  In  the  immediate  neighborhood  of  the 


Figure  30.  Variation  of  mean  Nus¬ 
selt  numbers  with  ambient  tempera¬ 
ture  Ta  (0,  experimental  results). 
The  solid  curve  is  from  the  analytical 
work  of  Gebhart  and  Mollendorf  (1978), 
the  dashed  curve  is  the  prediction  with 
the  Boussinesq  approximation  (after 
Bendell  and  Gebhart  1976). 
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flow  direction  inversion  7]v,  very  slow  flows  exist  with  the  effective  Grashof  number  become 
very  small,  so  the  validity  of  the  simplest  boundary-layer  theory  becomes  questionable;  no 
theoretical  results  for  this  regime  were  found  in  the  literature.  In  general,  the  theoretical  re¬ 
sults  of  Gebhart  and  Mollendorf  after  multiplying  by  a  factor  of  (0. 102/0.303)1/4  (solid  curve 
in  Fig.  30)  compared  remarkably  well  with  those  reported  by  Vanier  and  tien  (see  solid  curve 
in  Fig.  29),  even  though  a  rather  elaborate  and  more  accurate  density-temperature  of  water 
was  claimed  to  be  used  in  Gebhart  and  Mollendorf’s  study. 


FORCED  CONVECTIVE  HEAT  TRANSFER  OVER  A  MELTING  SURFACE 


Yen  and  Tien  (1963)  were  the  first  to  investigate  the  problem  of  laminar  heat  transfer  over 
a  melting  plate.  With  the  assumption  of  a  linear  velocity  profile  (v„  =  cy)  and  with  coordi¬ 
nates  fixed  on  the  melting  surface,  the  energy  equation  is  given  as 


cy 


dT 

dx 


+ 


ar 

dy 


d!T 

dy 2 


(44) 


where  vy0  is  the  interfacial  velocity.  Equation  44  was  solved  by  iteration  process  with  the  first 
approximation  by  letting  vyo  =  0  and  X  =  (c/9ax)uiy.  The  first  approximate  solution  is 


T-  T 

_1 _ SSL 

T-Tm 


1exp(-6J)d4 

0 


(45) 


- 

where  6  is  a  dummy  variable  and  a0  is  defined  as  J  exp(-6’)er.  The  final  solution  can  be  writ 
ten  as  0 


T-T, 


T  -T 

*  1  m 


]  exp(-6’  +  ~^\d6 

SSL  =  0 _ “N  / 


(46) 


00  . 

where  aN  is  the  limiting  value  of  the  sequences  of  a„  j  [exp(-6J  +  /36)/an^]d  as  n  —  oo.  The 
melting  rate  is  given  by  0 


in  which 

^  ^7"»-7m) 

<j>  _  - - - . 

qL 

The  Nusselt  number  is 

Nu  -(rjfef  •  <48> 

The  ratio  of  the  Nusselt  number  with  melting  to  that  without  melting  and  by  neglecting  fac¬ 
tor  (c/9ax)u}  is 


=  Nu/Nu„  =  fl0/crN. 


(49) 


If  a  reasonable  estimation  of  c  =  t^/h  is  given  where  rw  is  the 
wall  shear  stress,  the  ratio  of  the  Nusselt  number  becomes 


Table  3.  Numerical 
values  of  as  as  a  func¬ 
tion  of  0  (after  Yen 
and  Tien  1963). 

St  =  -  St;  a.v 


Since  r*  is  proportional  to  the  velocity  gradient,  and  with  the  as¬ 
sumption  that  the  distribution  of  the  velocity  profile  is  compara¬ 
ble  to  that  of  the  temperature  profile,  eq  50  becomes 

<H0)  =  Nu/Nu0  =  (a0/flN)4'3.  (51) 

Figure  31  shows  the  effect  of  0(  =  -Si)  on  the  values  of  v/(d)  and 
&(0).  Table  3  shows  the  values  of  aN  as  a  unique  function  of  0. 
From  the  numerical  values  of  aN,  it  becomes  obvious  that  the 
surface  heat  flux  reduces  the  steepness  of  the  temperature  gradi¬ 
ent  and  consequently  decreases  the  heat  transfer  coefficient.  Fig¬ 
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ure  32  illustrates  the  effect  of  0  on  the  temperature  distribution.  It  clearly  indicates  that  at 


the  same  value  of  X,  the  dimensionless  temperature  T-  Tm/Tx-Tm  is  lowered  as  0  is  increased 
from  0  to  1 . 


The  assumption  v„  =  cy  simplified  Yen  and  Tien’s  (1963)  analysis,  but  it  is  known  that  for 
fluids  of  high  Prandtl  numbers  the  thermal  boundary  is  not  as  thick  as  the  velocity  boundary 
layer.  In  view  of  this,  Tien  and  Yen  (1964)  performed  an  additional  analysis  of  the  same 
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Figure  31.  Relationship  between  4, 
h.  and  0. 


Figure  32.  The  effect  of  0  on  temperature 
distribution. 


problem  but  proposed  the  portion  of  the  velocity  profile  within  the  thermal  boundary  layer 
as 


v,  =  c(x)y. 

The  energy  equation  is  solved  as 


T'(q)  = 


T-T, 


T  -T 

1  30  1  m 


m_  _  0 


*  exp/-»/J  +  ~W») 
o  '  "n/ 


in  which 


=  ti 


dr_  | 
dq  I  0 


and 


V  =  |9o(^  [c(jic)]1/2 


The  Nusselt  number  is 


Nu  = 
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The  ratio  of  the  Nusselt  number  is  now 
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(52) 


(53) 


(54) 


(55) 


(56) 


(57) 


The  first  part  on  the  right  side  of  eq  57  is  equal  to  ^(/3)  =  aN,  and  if  it  is  assumed  that  the  ef¬ 
fect  of  melting  on  the  temperature  gradient  is  comparable  to  the  effect  on  the  velocity  gradi¬ 
ent,  then  eq  57  is  equivalent  to  eq  51,  indicating  that  results  obtained  based  on  a  rather  re¬ 
stricted  assumption  are  applicable  to  more  realistic  cases. 

Pozvonkov  et  at.  (1970)  also  conducted  an  analytical  study  on  heat  transfer  at  a  melting 
flat  surface  under  the  conditions  of  forced  convection  and  laminar  boundary  layer.  In  their 
study,  the  boundary  layer  equations  were  solved  with  the  assumptions  of  constant  thermal- 
physical  properties  of  the  fluid,  independence  of  the  ratio  of  momentum  to  thermal  boun¬ 
dary  layer  thickness  (6/<5t)  in  the  direction  of  flow,  and  a  fourth-degree  polynomial  represen¬ 
tation  of  the  velocity  and  temperature  distribution.  The  ratio  of  the  local  Nusselt  number 
with  melting  to  that  without  melting  is  given  as 


in  which  k,  is  the  Kutateladze  number  (  =  -  1  /St),  a  is  expressed  in  terms  kf  by 


where  T/Tm  and  u/um  are  the  dimensionless  thermal  and  velocity  boundary-layer  distribu¬ 
tions.  The  value  of  6,**  is  for  the  case  when  there  is  no  melting.  To  evaluate  the  values  of  6** 
and  6**,  the  values  of  a,J*  ( 6 ••  =  6**  when  there  is  no  melting)  has  to  be  evaluated  from 


6“  - d”  m) 

where  rj  =  y/ht  or  rj  =  y/b.  The  value  of  6g*  can  be  found  directly  from  eq  61,  but  the  value 
of  6,£*  will  contain  the  unknown  boundary-layer  thickness  ratio  e0  =  6t0/60.  By  retaining  on¬ 
ly  the  first  power  in  t0  in  the  expression  6,5*.  the  first  approximation  of  e0  can  be  evaluated 
from 

6<Ao  =  vPr  (62) 


for  the  case  of  no  melting.  Equation  62  is  reduced  from  the  general  expression  of 


in  which  \0  =  a/6kf  and  (  =  5,/6.  By  repeated  application  of  eq  61  and  the  expression  for 
6,**  derived  from  eq  60,  a  final  value  of 

6<Ao  (or  f  0  =  Sl0/b0) 

can  be  found.  This  final  value  of  f0  's  used  as  a  first  approximation  to  evaluate  the  values  of 
5**  and  6**  from  eq  60  and  61.  Equation  63  is  then  used  to  compute  a  new  value  of  6/6,  (or 
1/f)  and  is  used  again  in  eq  60  and  61  to  compute  a  second  set  of  values  of  6,**  and  6**.  A 
new  value  of  6/6,  is  computed  from  eq  63  and  is  repeated  until  this  final  value  is  nearly  identi¬ 
cal  to  the  previously  calculated  one.  The  results  from  this  rather  complicated  boundary-layer 
integral  analysis  are  compared  in  Figure  33  with  the  much  simplified  work  by  Yen  and  Tien 
(1963). 

The  effect  of  melting  on  forced  convection  heat  between  a  melting  body  and  the  surround- 


ing  fluid  was  studied  quantitatively  from  the 
points  of  view  of  boundary  layer  theory, 
film  theory,  and  penetration  theory  by  Tien 
and  Yen  (1965). 

In  the  case  of  boundary  layer  theory,  exact 
solutions  of  the  equations  of  motion, 
energy,  and  continuity  were  obtained  and 
the  relation  between  the  ratio  of  dT  ( =  h/hQ) 
of  heat  transfer  coefficient  with  melting  to 
without  melting  was  derived  in  terms  of  the 
parameter  i3(=  -  St)  as 

dj  =(r*'(0,Pr,0))/fi  (65) 


Figure  33.  Comparison  of  Nux/Nux0  vs  0. 


where  T *  is  the  derivative  with  >?  of  the  dimensionless  temperature  identically  defined  in  eq  53 
but  with  t)  defined  as 


V  =  Lv/2)(vjvx)'/2. 

The  parameter  K  is  given  as 
K  =  (0/Pr)T'  ’  (0,Pr,K). 


Numerical  values  of  [K  Pr/T*’  (0,Pr,K)]  as  functions  of  0T  have  been  compiled  by  Stewart 
(1950). 

In  the  analysis  based  on  film  theory,  the  resistance  to  the  interphase  transport  phenome¬ 
non  is  assumed  to  be  confined  to  a  thin  layer  of  stagnant  film  immediately  adjacent  to  the 
solid  boundary.  The  model  is  assumed  to  be  one-dimensional  (perpendicular  to  the  surface) 
and  of  steady  state.  The  parameter  0T  is  found  to  be 


ln(l  +  0) 


0 


(66) 


In  the  penetration  model  analysis,  it  is  hypothesized  that  the  transport  process  is  affected 
by  small  eddies  of  the  fluid  in  a  turbulent  field  coming  into  contact  with  the  interfaces.  Since 
in  most  cases  the  duration  of  contact  is  rather  short,  the  eddy  may  be  treated  as  a  semi-infin¬ 
ite  solid,  and  a  transient-state  heat  conduction  equation  can  be  used  to  describe  the  energy 
transfer  process.  A  solution  was  attained,  and  the  ratio  0T  is  expressed  as 


dj  =  //(l  —  erf  d>)  erf(d> ') 
in  which  <i>  is  related  to  0  by 


(67) 


(1 -erfO)  v7r«i>exp(<iJ)  =  ~0- 


(68) 


Figure  34  shows  the  effect  of  0  on  the  values  of  0T  from  the  considerations  of  boundary, 
film,  and  penetration  theory  along  with  the  works  reported  by  Yen  and  Tien  (1963)  and 
Merk  (1954).  Although  the  quantitative  results  differ,  they  all  show  the  same  qualitative 
trend  and  indicate  that  the  process  of  phase  transition  (melting)  inhibits  the  heat  transfer 
rate.  The  results  based  on  the  Leveque  solution  agree  with  those  based  on  boundary  layer 
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theory  for  small  values  of  $  but  deviate  from 
each  other  as  d  increases.  This  is  expected 
since  the  solution  obtained  by  Yen  and  Tien, 
in  a  sense,  is  an  asymptotic  solution  of  that 
based  on  boundary  layer  theory  for  small  ,i. 
The  difference  in  results  based  on  the  boun¬ 
dary,  film,  and  penetration  models  requires 
that  caution  be  exercised  in  selecting  among 
these  results  for  practical  applications. 

DISCUSSIONS  AND  CONCLUSIONS 

This  review  covers  only  the  problems  asso¬ 
ciated  with  the  anomalous  density-tempera¬ 
ture  relation  of  water.  The  discussion  and 
conclusion  are  grouped  into  three  sections, 
i.e.  1)  onset  of  convection,  2)  temperature 
structure  and  natural  convective  heat  trans¬ 
fer,  and  3)  laminar  forced  convective  heat 
transfer. 


Onset  of  convection  10  c  ^  i!j 

The  criterion  of  the  onset  of  convection  of 

a  water  layer  that  contains  a  density  extre-  Figure  34.  Comparison  of  various  theor- 
mum  was  found  both  experimentally  (Yen  'es:  r^e  effect  °f  &  on  (after  Tien  and 

1968,  Yen  and  Galea  1969)  and  analytically  Yen  1965). 

(Veronis  1963,  Sun  et  al.  1969,  Merker  et  al. 

1979)  to  be  not  a  constant  value,  as  in  the 

classical  B6nard  problem,  but  dependent  on  the  thermal  boundary  conditions  (if  the  layer  is 
formed  by  phase  transition,  then  one  of  the  thermal  boundaries  is  at  the  ice  melting  point, 
i.e.  Tm  =  0°C).  The  experimental  values  were  compared  with  analytical  values  very  favor¬ 
ably  as  shown  in  Figure  12  in  terms  of  Rac  experimental  vs  Rac  theoretical  and  in  Figure  1 3  as 
Rac  vs  7,  or  T2  explicitly  and  from  which  it  can  be  concluded  that  in  the  case  of  melting  from 
the  top,  the  higher  the  values  of  T},  the  greater  Rac  becomes,  or  in  other  words  the  farther  re¬ 
moved  temperature  7j  is  from  4°C,  the  less  prone  the  layer  is  to  the  onset  of  convection.  On 
the  other  hand,  in  the  case  of  melting  from  below,  as  temperature  7",  increases,  Rac  reduces 
exponentially  and  approaches  the  value  =  1708  asymptotically  as  reported  in  the  classical 
B6nard  problem.  This  is  evident  from  the  fact  that  as  T  becomes  higher  and  higher,  the 
buoyancy  forces  created  by  temperature  difference  A T  =  7",  -  Tmax  possess  a  much  stronger 
influence  on  the  layer  stability  than  the  effect  produced  by  the  density  extremum  (i.e.  = 
4°C),  and  subsequently  the  continuously  forming  layer  behaves  like  a  normal  fluid  (i.e.  there 
is  a  monotonic  density-temperature  relationship)  as  in  the  Benard  problem.  On  the  other 
hand,  as  T,  decreases  and  approaches  the  temperature  of  the  density  extremum  (  =  4°C),  Rac 
increases  and  approaches  the  limiting  value  of  infinity.  This  is  also  expected,  since  if  T,  is 
maintained  at  =  4°C,  the  water  has  its  highest  density  at  the  lower  boundary  and  the  water 
layer  will  always  remain  stable. 

In  the  case  of  a  water  layer  formed  by  mt  .ing  ice  from  above,  the  trend  of  variation  of  Rac 
with  boundary  temperature  is  just  reversed.  The  higher  the  temperature  T2,  the  greater  Rac 
becomes.  This  can  be  explained  by  noting  that  if  T:  is  maintained  in  the  range  0  <  7j  s  4°C, 
the  entire  layer  is  unstable  because  the  higher  density  water  will  overlie  the  less  dense  water 
underneath,  and  will  consequently  result  in  lower  Rat  values.  If  7j  is  maintained  at  a  higher 
temperature  than  4°C,  only  a  fraction  of  the  layer  (  =  4 H/T2)  is  potentially  unstable,  so  the 
layer  is  less  prone  to  onset  of  convection.  It  seems  that  the  Rac  value  grows  greater  and 

30 


■ft’ 

■  If. 


WWW 


greater  as  the  effect  of  the  density  extremum  becomes  less  and  less  pronounced.  It  is  also  in¬ 
teresting  to  note  that  the  two  Rac  curves  intersect  at  exactly  7",  =  Ti  =  8°C,  which  clearly  in¬ 
dicates  that  under  this  particular  thermal  condition,  these  two  systems  are  identical  and  have 
a  unique  Rac  value  regardless  of  how  the  water  layer  was  formed. 

Temperature  structure  and  natural  convective  heat  transfer 

The  most  striking  phenomenon  of  the  temperature  distribution  either  in  the  constant  water 
layer  depth  (Townsend  1964,  Myrup  et  al.  1970)  or  in  the  continuously  growing  layer  from 
melting  ice  (Yen  1984)  is  the  formation  of  a  nearly  constant  temperature  region  that  eventu¬ 
ally  expands  to  two-thirds  of  the  entire  layer  depth.  The  only  significant  difference  observed 
is  that  the  temperature  in  the  constant  temperature  region  in  the  water  layer  formed  by  melt¬ 
ing  ice  from  below  depends  on  the  boundary  temperature  T ,  (see  Fig.  19c),  but  is  a  fixed  val¬ 
ue  (=  3.2°C)  independent  of  the  imposed  upper  boundary  temperature  7)  (Fig.  21c)  in  the 
water  layer  formed  by  melting  from  the  top.  This  is  identical  to  the  temperature  structure 
measured  in  the  rather  deep  and  constant  depth  water  layer  reported  by  Townsend  (1964) 
(Fig.  17a)  and  Myrup  et  al.  (1970). 

The  heat  flux  for  melting  from  the  top  is  found  to  be  a  weak  function  of  T,  and  can  be  ex¬ 
pressed  as  q  =  177  (Ti)0303  in  which  q  is  given  in  W/m\  However,  in  melting  from  below, 
the  heat  flux  is  found  to  be  strongly  dependent  on  T,  and  can  be  approximated  by  q  = 
-  1900  +  315(7',).  This  is  because  as  T,  is  maintained  at  a  higher  value  there  will  be  a  stronger 
current  of  mixing  in  the  constant  temperature  region,  thus  reducing  the  laminar  layer  thick¬ 
ness  on  the  lower  boundary  and  the  upper  water/ice  interface  and  increasing  the  heat  trans¬ 
fer  rate.  On  the  other  hand,  the  effect  of  7)  on  heat  transfer  is  much  weaker  because  there  is 
a  stable  layer  overlying  the  expanding  unstable  region  and,  in  addition,  the  buoyancy  force  is 
created  merely  by  a  temperature  difference  of  4°C  (i.e.  4°C-0°C  =  4°C).  The  works  of 
Townsend  (1964)  and  Adrian  (1975)  are  similar  to  the  work  of  melting  ice  from  the  top,  and 
they  reported  nearly  identical  heat  fluxes  of  approximately  340  W/mJ,  which  is  somewhat 
lower  than  those  obtained  from  the  melting  experiment.  The  discrepancy  can  probably  be  at¬ 
tributed  to  the  unsteady  nature  of  the  melting  experiment  (i.e.  the  layer  depth  is  initially  at 
zero  and  it  deepens  as  the  melting  progresses). 

The  works  of  Tkachov  (1953),  Merk  (1954),  Schechter  and  Isbin  (1958),  Vanier  and  Tien 
(1968,  1970),  Bendell  and  Gebhart  (1976),  and  Gebhart  and  Mollendorf  (1978)  were  conduct¬ 
ed  either  theoretically  or  experimentally  and  were  aimed  at  classifying  the  implication  of  the 
density  extremum  on  the  heat  transfer  characteristics  of  a  melting  as  well  as  a  nonmelting 
system.  Tkachov  was  the  first  to  suggest  that  under  certain  thermal  conditions  the  boundary 
layer  might  be  split,  with  a  predominantly  upward  motion  immediately  adjacent  to  the  ice 
surface  and  a  region  of  downward  motion  outside  this.  Based  on  his  experimental  work  on 
ice  spheres  and  vertical  and  horizontal  cylinders,  he  reported  that  the  heat  transfer  coeffi¬ 
cient  is  lowest  at  a  7^  of  about  5.5°C.  Merk  was  the  first  to  take  up  this  moving  boundary 
problem  analytically;  he  solved  it  with  an  approximation  method.  For  the  case  of  no  melting 
and  with  a  small  value  of  Tx,  he  reported  the  inversion  temperature  to  be  at  =  5.005°C. 
Based  on  a  rather  complete  calculation  (i.e.  without  the  limitation  on  Tx),  Merk  reported  the 
inversion  temperature  to  be  around  5.30°C  for  no  melting  and  5.31  °C  with  melting;  he  fur¬ 
ther  indicated  that  the  effect  of  melting  is  only  appreciable  for  Tm  >  7]v  and  may  be  neglect- 
ed  for  Tm  <  Tjv  (Fig.  24). 

The  analytical  work  of  Vanier  and  Tien  (1968)  also  reported  the  existence  of  three  types  of 
flow  regions,  i.e.  laminar  upward,  downward,  and  dual  flow  (sometimes  termed  “inverted” 
or  “split”),  as  well  as  regions  of  no  solution.  These  regions  are  determined  by  specific  com¬ 
binations  of  rw  and  Tm.  Their  computed  results  are  found  to  be  consistent  with  the  observa¬ 
tions  by  Tkachev  (1953),  the  analytical  and  experimental  works  of  Schecter  and  Isbin  (1958), 


and  the  analytical  results  of  Merk  (1954)  as  well  as  the  experimental  results  of  Dumor^  et  al. 
(1953).  These  analytical  findings  were  verified  with  their  experimental  work  on  melting 
spheres  similar  to  those  reported  by  Tkachev  and  Dumor6  et  al.  A  minimum  heat  flux  oc¬ 
curred  around  7jv  =  5.35  ±0.2°C  as  compared  with  Merk’s  5.31  °C  and  Tkachev’s  5.5°C. 
Their  findings  are  also  in  good  agreement  with  the  most  recent  work  of  Bendell  and  Gebhart 
(1976).  Based  on  the  melting  of  vertical  ice  plates,  they  reported  an  inversion  temperature  of 
5.6°C.  Gebhart  and  Mollendorf  (1978),  however,  claimed  to  have  used  a  more  elaborate  and 
accurate  density-temperature  representation,  and  their  results  are  comparable  with  those  re¬ 
ported  by  Vanier  and  Tien  (1968),  even  though  a  less  complex  density-temperature  relation 
was  used. 

Laminar  forced  convective  heat  transfer 

Based  on  the  limited  analytical  work  reported  on  forced  convective  heat  transfer  over  a 
melting  surface,  it  can  be  concluded  that  the  interfacial  velocity  resulting  from  phase  transi¬ 
tion  tends  to  retard  heat  transfer.  From  the  most  simplistic  approach  to  this  problem,  Yen 
and  Tien  (1963)  found  the  ratio  of  the  Nusselt  number  with  melting  to  that  without  melting 
to  be  a  strong  function  of  0[=  cp(Tx-  T^/L  =  -St)  (Fig.  31).  This  result  was  in  good 
agreement  with  those  reported  by  Pozvonkov  et  al.  (1970)  from  a  rather  complicated  integral 
solution  of  the  transport  equations  (Fig.  33).  The  reduction  of  the  heat  transfer  rate  is  attrib¬ 
uted  to  the  lowering  of  the  temperature  gradient  and  its  equivalent  to  the  case  of  mass  injec¬ 
tion  through  the  laminar  layer  into  the  mainstream.  This  evidence  (Fig.  34)  was  further  dem¬ 
onstrated  by  Tien  and  Yen  (1965),  based  on  their  quantitative  analysis  of  this  moving  boun¬ 
dary  problem  with  classical  boundary,  film,  and  penetration  theory. 
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